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ABSTRACT 

We show that the cool gas masses of galactic discs reach a steady state that lasts 
many Gyr after their last major merger in cosmological hydrodynamic simulations. 
The mass of disc gas, Mgas, depends upon a galaxy halo’s spin and virial mass, but 
not upon stellar feedback. Halos with low spin have high star formation efficiency and 
lower disc gas mass. Similarly, lower stellar feedback leads to more star formation so 
the gas mass ends up nearly the same irregardless of stellar feedback strength. Even 
considering spin, the Mgas relation with halo mass, M 200 only shows a factor of 3 
scatter. The Mgas-M 2 oo relation show a break at M 200 = 2 x that corresponds 

to an observed break in the Mgag-M* relation. The constant disc mass stems from a 
shared halo gas density profile in all the simulated galaxies. In their outer regions, 
the profiles are isothermal. Where the profile rises above n = 10“^ cm“^, the gas 
readily cools and the profile steepens. Inside the disc, rotation supports gas with a 
flatter density profile and at low masses, supernova explosions disrupt the disc. Energy 
injection from stellar feedback also provides pressure support to the halo gas to prevent 
runaway cooling flows. The resulting constant gas mass makes simpler models for 
galaxy formation possible, either using a “bathtub” model for star formation rates or 
when modeling chemical evolution. 

Key -words: galaxies: formation - galaxies:ISM - hydrodynamics - methods: N-body 
simulation 


1 INTRODUCTION 

How gas evolves, especially gas at temperatures well below a 
galaxy halo’s virial temperature, remains one of the central 
questions of galaxy formation. The accretion and cooling of 
gas onto discs makes star formation possible. There are large 


reservoirs of hot gas surrounding galaxies (Tumlinson et al. 
|2011 Werk et al.|2014 l and centrifugally supported discs of 
cool gas that form stars. 

There are a number of physical processes that can move 
gas between the disc and halo. Processes that add gas to 
discs include: 

• Cold flow accretion of pristine cold gas through cold, 
filamentary accretion ( Keres et al.|2005 l 

• Hot mode accretion, which happens when gas radia- 


1977 [White fc Rees|p^978| [Putman et al.|[2012| [Werk et al.[ 


20141 


• Evolved stars: Stellar populations lose ~ 40% of their 
initial mass through stellar winds and planetary nebulae 
i Kalirai et al.|20^ Leitner fc Kravtsov||2011 1 

• Stimulated accretion: In a process not resolved in cur¬ 
rent cosmological simulations, metal enriched outflows from 
the disc can seed the halo with enough metals and turbu¬ 
lence to encourage accretion onto the disc i jFraternali fc Bin^ 
ney||2008 Marinacci et al.|2010 | 

• Gas recycling: In a similar process to stimulated accre¬ 
tion, simulations show that metal rich gas ejected through 
stellar winds often falls back onto the disc ||Oppenheimer 
et al.|2010 Brook et al.|2014 Ubler et al.|20lj | 


tively cools out of the hot halo onto the disc ( [Rees fc Ostriker| Sinks that can remove gas from discs include: 


* Email: stinson ‘at’ mpia.de 


• Star formation 

• Galactic outflows: stellar or AGN driven 
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2 Stinson et al. 


• Photoionization or cosmic rays: strong photo ionizing 
sources photo-evaporate cool gas (Kannan et al.|2014l 


With such a wide variety of processes that add and re¬ 
move cold gas from galactic discs, it seems like the disc gas 
mass would not stay constant. However, there have been a 
number of recent steady state models of galaxy formation 


that suppose exactly that (Bouche et al.|2010 

Bauermeister 

et al.|2010[ 

Dave et al.|2012[ [Lilly et al.|2013| 

Dekel & Man- 

delker|2014 

). Such models present a relatively simple picture 


of galaxy formation and with such a simple description, they 
do a remarkably good job reproducing key observables like 
the evolution of the specihc star formation rate. It is thus 
worthwhile exploring in numerical simulations: Does the as¬ 
sumption of constant gas mass hold? 

The dynamical time for the inner regions of galaxies 
immediately surrounding the disc is only a few percent of 
the Hubble time. So, in the absence of mergers, it seems 
plausible that all the processes listed above might balance 
to create a quasi-steady state. 

There are a variety of observational surveys of the cold 
gas content in galaxies that provide a basis for comparison 
with models. For neutral hydrogen (Hi), there have been the 
HI Parkes All-Sky Survey hipass ||Meyer et al.||2004|), the 


Arecibo Legacy Fast ALFA (alfalfa) Survey (Giovanelli 


et al. 20051, and the Galex Arecibo SDSS Survey (gass) 
( [Catinella et al.]|2012[ ). The COLD GASS Survey augmented 
GASS with the inclusion of CO observations that give some 
idea about the mass of hydrogen in the molecular phase. 
These surveys of the local Universe provide a good census 
of the cold gas in galactic discs. 

There have been observations of CO throughout cosmic 
history back to z ~ 3 (Daddi et al.|2010 Genzel et al.|2010 


Tacconi et al.|[^13). However, the 21 cm emission from Hi 


is too faint to see beyond the local Universe. For this rea¬ 
son, more powerful instruments like the Square Kilometer 
Array (SKA) ( Blyth et al.||20f5 l are required to nail down 
the total gas mass evolution of galaxies. Using the spectra 
of absorption lines from distant quasars. [Prochaska et alT| 
(20051 found that the cosmic density of Hi, flui, is constant 
(also see |Crighton et S(]|2015[ ). 

Relatively simple models have been able to match the 
data well. Semi-empirical methods estimate the gas con¬ 
tent of galaxies based on derived star formation rates from 
abundance matching ( Lu et ar]|2015 Popping et al.pOlS l. 
Such a model extrapolates H2 masses from the star forma¬ 
tion rate based on the Kennicutt-Schmidt star formation- 


gas surface density relationship (Schmidt 1959 Kennicutt 
1998 Bigiel et al.|2008 L The model further extrapolates Hi 


masses from gas pressure arguments (Blitz & Rosolowsky 


2006). The model does a good job matching H2 masses at 
redshifts z ~ 2. These analytic models imply little evolu¬ 
tion of the relationship between Mhi and Mvir across cosmic 
time, though the exact location of that gas unknown. Some 
may be in discs, but much could be in the cosmic web or 
surrounding galaxies. 

Similar estimates have been applied to other galaxy for¬ 
mation models to show the gas content and evolution of a 
wide range of galaxy models. Using a model for galaxy for¬ 
mation based on angular momentum conservation, |Dutton| 
|et al.| ( |2010[ ) found constant cold gas mass with a grad¬ 
ual shift from molecular to atomic hydrogen as galaxies 


grow larger and their surface densities decrease. [Lagos et al.| 
(20111 showed how gas evolves using the GALFORM semi- 


analytic model from Durham. Popping et al. (2014) exam¬ 


ined gas properties of the Somerville & Primack 


(19991 semi- 


analytic model and found a similar result. [Popping et al.j 


(2015) and Lu et al. (2015) used a semi-empirical approach 


to constrain the gas evolution in galaxies much like has been 
done with stars using abundance matching and halo occu¬ 
pation distributions. In their semi-empirical work, [Popping 
et al. (2015) points out that the Hi content of galaxies does 


not evolve much, exactly like we emphasize here. 

Numerical simulations of galaxy formation have entered 
an era in which it is possible to create large samples of re¬ 
alistic galaxies (Aumer et al.|20T3 Hopkins et al.|2014 Vo- 


gelsberger et al. 2014 Schaye et al. 2015 Wang et al. 1 201 5' ). 

Genel et al.| ( |2014p uses a simple formulation to estimate 
the H2 fraction in the Illustris galaxies and hnds agreement 
with observations at 2: = 0, but not at 2 = 2. Some studies 
have focused on the gas content and distribution in galaxy 
simulations. In a paper on the stellar mass metallicity rela¬ 


tion, Finlator & Dave (2008) found minimal changes in disc 
gas mass and [Dekel et al.| ( 2013P found the same in a pa¬ 
per focused on comparing simulations with analytic models. 
Nelson et al. (2015) examines several zoom regions at 2 ~ 2 


and finds that gas distributions are asymetric, particularly 
as resolution increases. 

We examine the cold gas properties and evolution in 
Numerical Investigation of a Hundred Astrophysical Ob¬ 


jects, NIHAO (Wang et al. 2015). The NIHAO project 


strives to simulate 100 galaxies across a wide mass range, 
10® <M*/M0 < 10®®, using hydrodynamics in a cosmologi¬ 
cal context. The NIHAO simulations are based on the star 
formation and feedback that was developed for the Mak¬ 
ing Galaxies in a Cosmological Context project (MaGICC, 
Stinson et al.|[2013 ). The MaGICC simulations used physi¬ 
cally motivated feedback to create a realistic galaxy of about 
Milky Way mass. The NIHAO project extends the range of 
simulations down to the smallest mass galaxy that will form 
stars, ~ 1O®M0. Using the same stellar physics across the 
entire mass range, [Wang et al.| ( |2015| ) find that the stel¬ 
lar mass of all galaxies agree with results from abundance 
matching. This makes the NIHAO sample ideal for exploring 
the physics of galaxy formation. 

Here, we examine the gas, especially the cold gas, in the 
ensemble of the NIHAO galaxies. We Hnd that despite all the 
complicated physics, the cool disc gas mass remains largely 
constant within a mass range of 10® <M*/M0 < 10®®'® and 
at moderate spin parameters, 0.03 < A < 0.06. We explore 
what the ramihcations of the constancy might be for study¬ 
ing galaxy evolution. In we provide an overview of the 
simulations used to find this result. shows the constant 
gas mass and how it relates to halo virial mass. !|^ discusses 
the gas density profiles that form in the simulations, ^pro¬ 
vides a context for this finding in recent theoretical models, 
observations of gas in galaxies, and some possible implica¬ 
tions for future observations. 


2 SIMULATIONS 

The simulations analyzed in this paper are taken from the 


NIHAO project (Wang et al. 2015). The NIHAO galax- 
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Constant Disc Gas Mass 3 


ies come from (20 h ^ Mpc)^ and (60 h ^ Mpc)® cubes 


from Dutton & Maccio ( 2014|). T hey use parameters from 
Planck Collaboration et al. ( |2014 ): Din=0.3175, flA=0.6825, 


f2bary=0.049, Ho = 67.1 kmT^pc■^ ag = 0.8344. 


2.1 ESF-Gasoline2 

In their paper describing superbubble feedback, [Keller et al.| 
| |2014[ ) described the updates they made to the N-body SPH 
solver GASOLINE (Wadsley et al.|2004l. While we do not use 


c* is the efficiency of star formation, i.e. the fraction of gas 
that will be converted into stars during tdyn- 

The simulations use the same feedback described in 
[Stinson et al.| ( [2013[ ). In this scheme, the stars feed en¬ 
ergy back in two epochs. The first epoch, “pre-SN feed¬ 
back” (ESP), happens before any supernovae explode. It 
represents stellar winds and photoionization from the bright 


young stars. In GASOLINE, as in Stinson et al.| (20131, the 


their superbubble feedback, we employ their modihed ver¬ 
sion of hydrodynamics that removes spurious numerical sur¬ 
face tension and improves multiphase mixing. We thus refer 
to our version of GASOLINE as ESF-GASOLINe2, where esf 
stands for the “early stellar feedback” that is described be¬ 
low. The biggest differences to the hydrodynamics in ESE- 


pre-SN feedback consists of 10% of the total stellar flux be¬ 
ing ejected from stars into surrounding gas (2 x 10®° erg of 
thermal energy per M© of the entire stellar population). We 
show ESF-GASOLINe2 simulations that use eBSF=10%. Be¬ 
cause of the increased mixing in ESF-GASOLINe2, the sim¬ 
ulations required more stellar feedback to have their star 
formation limited to the abundance matching value. Thus, 
we also show results from using eFSF=13%, which gives a 


better match to Behroozi et al. (20131 abundance matching 


2001 

1 proposed for calculating ^. 

Ritchie & Thomas 

2001 


also proposed modifying the density calculation to use equal 
pressures, but we do not use those densities in the simula¬ 
tions described here. 

Diffusion of quantities like metals and thermal energy 
between particles has been implemented as described in 


Wadsley et al. (2008L Metal diffusion is used, but ther- 


results. Radiative cooling is left on for the pre-SN feedback. 

The second epoch starts 4 Myr after the star forms, 
when the first supernovae start exploding. Only supernova 
energy is considered as feedback in this second epoch. Stars 
8 Mq < M* < 40 Mq eject both energy and metals into 
the interstellar medium gas surrounding the region where 
they formed. Supernova feedback is implemented using the 


mal diffusion is not used because it is incompatible with 
the blastwave feedback that delays cooling. ESF-gasoline2 
includes several other changes to the hydrodynamic calcu¬ 
lation. The Saitoh & Makino (20091 timestep limiter was 


blastwave formalism described in Stinson et al. (2006 L Since 


implemented so that cool particles behave correctly when a 
hot blastwave hits them. To avoid pair instabilities, ESF- 
GASOLINE2 uses the Wendland C2 function for its smooth¬ 
ing kernel ( Dehnen fc Aly|2012 I. The treatment of artificial 
viscosity has been modified to use the signal velocity as in 


Price (20081. 


the gas receiving the energy is dense, it would quickly be 
radiated away due to its efficient cooling. For this reason, 
cooling is delayed for particles inside the blast region for 
~ 30 Myr. 


2.3 Cool Gas versus HI+H 2 

We use two different measures of the cool gas content in the 
paper, one for comparison with observations {MHI+H 2 ), and 


The cooling is as described in IShen et al.l poTo| and is more physically relevant, M3,,('^20000 K). 


was calculated using CLOUDY (version 07.02; Ferland et al. 
([1998P) including photoionization and heating from the 


Haardt & Madau (20051 UV background, Compton cool- 


For comparison with observations (Fig.[^, we show the 
non-equilibrium calculated Hi.The non-equilibrium cooling 
calculation in the simulations tracks the evolution of the 


ing, and hydrogen, helium and metal cooling from 10 to 10° 
K. In the dense, interstellar medium gas, we do not impose 
any shielding from the extragalactic UV field as the extra- 
galactic field is a reasonable approximation in the interstellar 
medium. 


2.2 Star Formation and Feedback 

The simulations use a common recipe for star formation de¬ 
scribed in Stinson et al. ( [2006 1 that we summarize here. 
Stars form from cool {T < 15,000 K), dense gas. The metal 
cooling readily produces dense gas, so the star formation 
density threshold is set to the maximum density at which 
gravitational instabilities can be resolved, {rith > 9-3 

cm 


atomic ion species of hydrogen and helium (Shen et al. 20101. 


We report the neutral hydrogen mass as H 1 -I-H 2 . The simu¬ 
lations do not resolve densities high enough for H 2 to form, 
so we are not able to subtract the H 2 mass from our neutral 


gas masses. Saintonge et al. (20111 shows that Mh 2 are at 


where Mgas = 2.2 x 10® Mq and e is the gravitational 
softening (310 pc). Such gas is converted to stars according 
to the equation 


AM* 

At 




c* 


tdy 


( 1 ) 


Here, AM* is the mass of the star particle formed. At is the 
timestep between star formation events, 8 x 10® yr in these 
simulations, and tdyn is the gas particle’s dynamical time. 


most equal to Mhi, so in some cases, it is possible that our 
masses should be divided by a factor of 2. 

Hi observations only represent emission from the 21 cm 
transition that only exists at column densities above 10^® 
cm“^. It is possible that there is a small contribution to the 
simulated Mhi from lower column density gas, but the UV 
background should ionize that gas. 

As a measure of the total virialized halo, we use every¬ 
thing contained within the surface where the mean density 
is 200pc. We refer to this radius as rooo and the mass inside 
it as M 200 . We And that M(T< 20000 K) is more tightly cor¬ 
related with M 200 than Mhi-i-H 2 , so for all plots other than 
Fig. a we report our findings in terms of M(T<20000 K). 


3 RESULTS 

Our study of the gas mass in galaxies starts with a compar¬ 
ison with observations to check whether a detailed analysis 
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Figure 1. Simulated Hi masses as a function of stellar mass com¬ 
pared with an observed sample of galaxies from [Radford et al.| 
j |2015[ |. The simulated points are coloured according to the total 
spin parameter at 2 = 0 of the entire galactic halo. Low spin (red) 
galaxies have lower gas to star ratios than higher spin galaxies be¬ 
cause they are more efficient at turning gas into stars 


is worthwhile. Since the correspondence with observations is 
good, we investigate how the gas mass evolves and examine 
the relationship with halo mass. We also try to inform the 
trend we sees with some tests of galaxies as a function of 
spin parameter and using various recipes for feedback. 


3.1 Comparison with Observations 

Fig .[^ shows how the present-epoch masses of cold (H 1 -I-H 2 ) 
gas discs in the simulations compare with observations. The 
simulation data are shown as the points coloured according 
to spin. The grey points are observations from the|Bradford| 


et al. (20151 sample of isolated galaxies. The grey line shows 


the fit Bradford et al. (20151 make to their observations. The 
fit is a double power law with a break at M* = 10® ®Mq. Be¬ 
low the break, Bradford et al. (20151 find Mhi+H 2 c< 
above the break, they find Mhi+h 2 oc 

The [Bradford et al.| ( |2015| ) Hi masses come from 21- 
cm emission lines of their own low mass galaxy sample and 
their re-measnrements of the alfalfa survey 40% data re¬ 
lease ( [Haynes et al.|2011| ). These Hi observations have been 
cross matched to a re-rednction of the Sloan Digital Sky 
Survey (SDSS DR8), where stellar masses have been mea¬ 
sured with kcorrect using SDSS ugriz photometry, GALEX 


UV photometry and a Chabrier (20031 IMF. We have also 


included their Mh 2 mass estimates, which are derived from 
Ha star formation rates. See [Bradford et al.[ ( [2015[ ) for a 
complete description of these three mass estimates and the 
isolation criteria they applied to this sample. 

The cold gas masses from the simulation are derived us¬ 
ing the non-equilibrium cooling calculation performed in the 
simulations. The correspondence between the observations 
and simulations is notable, especially the break in the rela¬ 
tionship. Both the simulated and observed relation steepen 


near M* < 10® Mq (a break first noted in Huang et al. 
[2012[ ). In general, the scatter of the observations is larger 
than the simulations. 

There are observed galaxies with higher Hi masses than 


the simulations at all stellar masses. Both the lBradford et J7I 
(20151 sample and the alfalfa survey are flux limited and 
therefore biased in detecting the most gas-rich galaxies in 
the local universe. At low M*, the simulations exhibit sys¬ 
tematically lower Mhi+H 2 than the observations. The obser¬ 
vations become less complete at these lower masses. 

Another notable trend in Fig.[^is that the low spin (red 
points) galaxies tend more towards the lower-right of the 
observations. The lower right on this plot represents higher 
stellar mass and lower gas mass, thus the gas to star ra¬ 
tio is lower for these galaxies. In other words, gas is more 
effectively converted to stars if the galaxy has a low spin 
parameter (also see [Huang et al.|2012[ ). |3.3[ expands on this 
finding. 


Catinella et al. (20121 point out that around 30% of 


galaxies at all stellar masses have no detectable 21 cm emis¬ 
sion from Hi in their Hi blind GASS Survey. Many of these 
“gas free” galaxieL] are in groups or clusters where their 
Hi has been stripped. Our simulated selection is entirely 
composed of isolated galaxies that will not experience these 
environmental effects. Thus, every high mass galaxy in the 
simulated sample contains observable gas, similar to what 


Bradford et al. (20151 observe in their isolated sample. The 


selection is necessary for simulations since simulating the 
hydrodynamics of galaxies flying through hot cluster gas is 
too computationally demanding at present to maintain the 
resolution of the individual galaxies. 

Fig.0 shows that the correspondence between simula¬ 
tions and observations is encouraging enough to lend cre¬ 
dence that the simulations can be further analysed to see 
how gas evolves in real galaxies. 


3.2 Disc Gas Mass Evolution 

Fig. shows the evolution of the cool gas (T< 20000 K, 
green line, not the same as Mhi shown in Fig.[^ mass inside 
Vgai = 0.2r2oo for a sample of ten NIHAO galaxies. The left 
column shows galaxies for which the disc gas mass remains 
within 20% of the median since the end of the era of major 
mergers at 2 ~ 1.5. The right column shows galaxies in 
which the gas mass evolves gradually, though one should 
notice that the evolution is less than 50% from the median. 
All of the gas mass evolution plots have been normalized 
so that their median value over this period of evolution is 1. 
The range of the y-axis only spans a factor of 4 and all of the 
galaxies that evolve stay within that range. The sample of 
galaxies that maintain a constant gas mass are significantly 
constant. Their variation from their median values is never 
more than 20%. 

To see how the variation in gas mass compares with 
the variation of the total halo mass, the total halo mass is 
shown as the blue line. The total halo mass has also been 
normalized so that its median is 1. The total mass is gen¬ 
erally 2 orders of magnitude higher than the gas mass. The 
evolution of the total gas mass shows when mergers occur. 
The bottom two constant cases show non-negligible merger 


^ The quotes around “gas free” refer to the fact that while there 
is no observable dense cool gas, but there is plentiful lower density 
hot and cold gas in the galaxy halos surrounding the stars (i.e. 
Werk et al.|2014 1. 
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Figure 2. Disc gas mass evolution for ten high resolution galaxies from project NIHAO [green line). The values have been normalized 
by their median value. The dashed cyan line is a linear fit to the evolution of the gas mass from the time of its last major merger until 
2 = 0. For comparison, the evolution of the total halo mass is shown as the solid blue line. The galaxies on the left are a representative 
sample across the mass range of galaxies that maintain their constant disc mass, while the right column shows galaxies in which the disc 
gas mass gradually evolves. 


histories. The cool gas mass evolves less than the total mass 
in both of these cases. 

The gas and virial masses have been included along 
with the spin parameter, A, at 2 = 0 for each of the 
galaxies. The spin parame ter, A. shown is de fined using 
A' = J/\J‘^GM^QQr 2 oo from Bullock et al. (20011 for all the 
particles, both dark and baryonic, inside r 2 oo- Both selec¬ 
tions contain galaxies from nearly the entire mass range, 
though the evolving selection (right panels) includes higher 
mass galaxies. It is worth exploring what sort of galaxies 
maintain constant gas masses and which ones evolve. 

To quantify the evolution of gas masses, linear fits to the 


gas mass evolution are made from their last major merger 
until 2 = 0. A major merger is defined as the dark matter 
mass changing by more than 30% over 3 snapshots. The fits 
are shown in Fig. as dashed cyan lines. 

The slopes of the linear fits for the entire NIHAO sample 
are shown in Fig. The slopes are plotted as a function of 
halo mass (left) and spin parameter (right). In the halo mass 
panel, the point colour represents the halo’s spin parameter 
at 2 = 0. In the spin parameter panel, the point colour 
represents the halo mass. 

The grey region shows the region in which gas mass 
changes by less than order unity over the Hubble time. 
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Figure 3. Disc gas mass evolution from every high resolution galaxies in project NIHAO as a function of halo mass (left panel) and 
2 = 0 spin parameter (right panel). Mgas represents a linear fit to the evolution of cool disc gas mass between 2 = 1 and 2 = 0. Those 
slopes are normalized by the median gas mass during this period to put galaxies of varying masses on even footing. 


Most points lie comfortably with this region, emphasizing 
the point that gas masses remain constant for significant 
time periods. 

Both panels have scatter. The points represent galaxies 
with a wide variety of merger histories. As we will see, gas 
mass typically follows total halo mass, so mergers tend to 
increase the disc gas mass of galaxies. Late mergers can also 
change the spin parameter at the last moment, so A mea¬ 
sured at 2 = 0 may not be the best time to measure the halo 
spin, but it is the one that could most readily be observed. 


Despite the complex formation histories for the galax¬ 
ies, it is still possible to pick out a trend with halo spin: the 
only galaxies that lose disc mass have Az=o < 0.042. That 
means galaxies that spin less lose cool gas content faster than 
those that are spinning quickly. One can understand such a 
reduction in cool gas with a simple model. Cool gas is ro- 
tationally supported. In galaxies with less spin, there is less 
rotational support, so the cool gas collapses towards the cen¬ 
tre of the galaxy disc. As it collapses, it reaches higher sur¬ 
face densities. Gas turns into stars more readily as it reaches 


higher surface density (Schmidt 1959 Kennicutt 1998 Bigiel 
et al.|[2008 1. Thus, the cold gas mass starts to decrease. It 


is notable that the cosmic median spin parameter is ~ 0.04 
( Bullock et al.|200T Maccio et al.|2008 l, the spin parameter 
below which galaxies decrease in cold gas mass. 


There is no trend apparent in cool gas mass as a func¬ 
tion of halo mass (left panel). At all halo masses, there is 
roughly an equal number of galaxies that lose cool gas mass 
to those that gain it. We might expect a decrease in cool gas 
masses at M200 > lO^^M© because the virial temperature 
rises high enough that it takes longer for gas to cool. Such a 
decrease would suggest how star formation can be quenched. 
Our sample of galaxies does not probe high enough masses to 


make such a statement. In the sample of galaxies at IO^^Mq 
there is again an even distribution of galaxies gaining and 
losing mass. One galaxy at 4 x has decreasing gas 

mass but does not make a trend. 


3.3 The Effect of Spin Parameter 

Following up on the spin trend from Fig.[^ we make a closer 
examination of how spin can affect disc gas mass. We turn 
to the simulations from Herpich et al.| ( |2015[ ) that are inde¬ 
pendent of galaxy merger history. In their study of isolated 


collapsing spheres of gas, Herpich et al. (20151 show the halo 


spin parameter sets the surface density profile for galactic 
discs. Fig. shows another result of [Herpich et ah] | |2015[ ), 
namely, the gas fraction of the galaxies depends on how fast 
the initial sphere is spinning. The dot-dash lines show the 
evolution of the total baryonic (gas and stars) mass while 
the solid lines show the evolution of the gas mass alone. The 
range of spin parameters shown (0.02 < A < 0.06) spans 
the middle of the range of spin parameters found in simu¬ 
lations of cosmological volumes (Bullock et al. ][200l| . Even 
in this small range, one can see the effect spin has on disc 
gas evolution: First, the variation in spin parameters leads 
to a factor of 3 spread in final gas mass. Second, the high 
spin galaxies continue to add gas throughout their evolution 
while the disc gas mass decreases in the low spin galaxies. 
Such results are consistent with the analysis injMo et al.j 


(19981 that shows how the scale radius of discs depends on 


the spin parameter and virial radius, 


Rd oc Ar 2 oo- 


( 2 ) 


The evolution of the galaxy sizes is the result of the ac¬ 
cretion of higher angular momentum material. Pristine gas 
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Figure 4. The solid lines show the evolution of the baryonic disc 
mass in isolated (non-cosmological) simulations of 10^^ Mq halos 
that include spinning gas spheres that cool and collapse to form 
stars. The dashed lines show the total baryon content in a disc 
region with a radius of 20 kpc and 6 kpc thickness. Because all 
the halos have the same mass, it is not surprising to see that the 
same amount of material cools into the disc region. There is an 
apparent difference in how much of that material remains as gas. 
The gas that is given the most initial spin (A = 0.06) supports the 
most gas and the halo with the least spin (A = 0.02) supports the 
least gas. The high spin galaxies also increase in gas mass while 
the low spin galaxies decrease in cool gas mass. 
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Figure 5. The cool disc gas mass as a function of M 200 at 4 
different epochs, 2 ~ 4, 2 ~ 2, 2 ~ 1, and 2 ~ 0. The same 
two lines are drawn on each panel. At low masses, the dotted 
line represents the relationship Mgas oc -^200’ while at M 200 > 
2xlO^^M0, Mgas = 0.04M200- The points are coloured according 
to the spin parameter of the galaxy at the moment the mass is 
measured. 


At M 200 > 2 X IO^^Mq, the relationship does not evolve 
much, but is consistently shallower than Mgas oc M 2 oo- 


that gets accreted onto the disc at late times has taken the 
longest time to get there because it started the furthest 
away. Thus, it has more angular mome ntum than t he gas 


that was accreted earlier. 


Brook et al. 


(2012a) and 


Ubler 


et al. (20141 point out that gas ejected from galaxies can 


also be re-accreted. Their studies show that ejected gas con¬ 
sists initially of the lowest angular momentum material and 
when it is re-accreted, it has more angular momentum. 


3.4 Full Population Statistics 

We have seen disc gas masses remain mostly constant over 
the course of many Gyr of evolution. There are some fluctu¬ 
ations and systematic trends depending on the spin parame¬ 
ter but not halo mass. To gain a more complete understand¬ 
ing of what sets the cool disc gas mass, Fig. shows the 
cool disc mass as a function of total halo mass at 4 differ¬ 
ent epochs. The points are coloured based on the halo spin 
parameter measured at the same time as the cool disc gas 
mass. At each epoch, there is a remarkably tight relationship 
between cool gas mass and virial mass that does not evolve 
much between epochs. No correlation with spin parameter 
is immediately apparent. 

To see how these relationships evolve, the same two lines 
are drawn in each epoch panel. At M 200 < 2 x IO^^Mq, 
the dotted line represents Mgas oc M 200 , while at M 200 > 
2 X IO^^Mq, the dashed line shows Mgas = 0.04M2oo- For 
M 200 < 2x IO^^Mq, the relationship between M 200 and Mgas 
evolves with time; it becomes steeper as the galaxies evolve. 


3.5 The Effect of Hydrodynamics 

While Fig. looks at a sample of galaxies simulated with 
uniform physics, it is interesting to examine what happens 
when the physics changes in a limited sample of two simu¬ 
lated galaxies. For this study, we use two cosmological zoom 
simulations drawn from the McMaster Unbiased Galaxy 


Simulations (MUGS, Stinson et al. 20101. The MUGS ini¬ 


tial conditions were created using a WMAP3 cosmology 
(f2m=0.24, Ua=0.76, flbary=0.044:, Ho = 73 km s“^Mpc“^, 
as = 0.76 Spergel et al.|20(i7 1. See Stinson et al. (20101 for a 
complete description of the creation of the initial conditions. 

To study which aspects of the numerical implementa¬ 
tion affect the disc gas mass. Fig. shows how the disc gas 
mass evolves when stellar feedback is varied or when hydro¬ 
dynamics changes. The definition for cool gas mass is slightly 
different in this case. The gas mass is measured inside a disc 
region 20 kpc in radius and 6 kpc in height centred on the 
disc midplane (the disc region includes 3 kpc above and 3 
kpc below the midplane) rather than a ~ 40 kpc sphere. In 
Fig. [I the solid lines represent the total baryon mass and 
the dashed lines represent only the gas portion of that mass. 

The two galaxies studied here have masses of ~ 7 x 10^^ 
M 0 , but different merger histories. In both cases, the discs 
reach a maximum gas mass and remain relatively constant 
for the rest of the simulation. In gl536, the change of hydro- 
dynamic scheme from GASOLINE to ESF-GASOLINe2 causes 
a change in the disc gas mass, but a change in feedback en¬ 
ergy in ESF-GASOLINe2 causes no change in disc gas mass 
even though the mass of stars nearly doubles in the lower 
feedback case. While we show a limited sample here, hy- 
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Figure 6. The Disc Gas Mass Conspiracy: A comparison of how the baryonic disc mass evolves using different hydrodynamic models 
for the galaxies. The solid line represents all the baroynic mass while the dashed line represents only the gas. Note that the y-axis is a 
linear axis in units of 10^*^ Mq. Cool gas mass is mostly independent of hydrodynamic and feedback treatment. In all cases, the cool gas 
mass reaches a nearly constant value. 


drodynamics likely plays a significant role in determining 
how much gas mass a disc can support. Thus, we propose 
that future studies might use disc gas mass diagnostic for 
comparing various hydrodynamic schemes. 

The largest difference in feedback is shown in the MUGS 
version of the simulations. In the MUGS simulations, the 
feedback was much weaker because MUGS employed a less 
top-heavy IMF and used only 40% of the supernova energy 
for feedback. Thus, stars output three times less energy per 
stellar mass in MUGS than the MaGIGC simulations. More 
than twice as many stars form in MUGS and the MUGS 
simulations contain less cold gas than any of the others. 
In the case of gl536, though MUGS forms so many more 
stars, the mass of cool gas is only off by ~ 50%, and the gas 
mass stays relatively constant throughout the simulation. 
In g5664, the gas mass is more like 30% the gas mass of 
the simulations with the other physics, but it again remains 
largely constant. 


4 UNIVERSAL HALO GAS PROFILE 

Fig. [5] shows that the cool gas mass is well correlated with 
the halo virial mass. Fig. shows that a range of feedback 
treatments on the same initial conditions result in more sim¬ 
ilar disc gas masses than stellar masses. With sufficient feed¬ 
back, the disc gas masses remain constant over many Gyrs. 
Gombining these facts indicates that the gaseous halos of 
galaxies might reach some state of hydrostatic equilibrium. 


4.1 Why is the disc gas mass regulated? 


To understand the steady state gas mass, it may be helpful 
to think of galaxy formation in the context of the hydro- 
dynamic balance between pressure and gravity (for a more 
complete description, see §4 in [White fc Fre^|1991[ ). The 
mass of material in a galaxy’s halo determines how much 
gravity is pulling gas into the center of the halo. The hydro- 
dynamic pressure provides support for the gas. The balance 
between these forces configures the gas into something like 
an isothermal sphere with a density profile of p oc r~^. As 
the density increases to small radii there is a radius called 
the cooling radius, r^, inside which the gas cools on a shorter 
timescale than the dynamical time. 


Field et al. (19691 described the instability of this con¬ 


figuration. Gas that cools removes its pressure support, so 
more gas flows to fill the vacuum, becomes dense and cools. 
What can limit this instability is some sort of pressure or 
energy support that keeps the halo gas at low enough densi¬ 
ties so that it does not cool too quickly. In our simulations, 
stellar feedback plays this role. The ultimate fuel source for 
stellar feedback is the nuclear fusion going on in stars. How¬ 
ever, to create the stellar feedback, stars must be formed, 
which requires cold gas, which leads to the steady state gas 
mass. 

To explore how that equilibrium looks. Fig. shows 
spherically averaged density profiles for all the NIHAO 
galaxies. Each point in the profile is coloured according to 
the median temperature of the gas in that radial bin. While 
not all the profiles look exactly the same, they share 3 com¬ 
mon sections: the disc, a transition, and an outer halo. To 
understand the profile shape, we examine each part sequen¬ 
tially, following the flow of gas from the outside in. 
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Figure 7. Density profiles for all the gas in the NIHAO sam¬ 
ple split into 4 different halo mass ranges. Each density point is 
coloured according to the median temperature for the radial bin. 
The radii at which the densities are measured are evenly space 
with the log of the radius. 



Figure 8. Gas density profile from r2oo iiito where the gas den¬ 
sity surpasses 10“^ cm~^. The profiles are coloured according to 
the the mass (M 200 ) of their halo. There are 100 logarithmically 
spaced radial bins. 


Fig-i focuses on the outer profile. The outer halo gas 
has a nearly isothermal profile with a slope close to -2. The 
outer profiles line up neatly according to their halo mass. 
Higher mass galaxies have higher virial radii, so the profiles 
have higher densities at radii more distant from the centre. 
All the galaxies in the sample follow a nearly isothermal 
shape up to a density of 10“® cm“®. 

After the radial profiles reach a density of 10“® cm“®, 
the profile transitions to a significantly steeper slope. Over 



log(n [cm ®]) 


Figure 9. Gas temperature-density phase diagram for gl.92el2, 
a M200 = 2.34 X 10 "^ Mq halo. The diagram is split into four 
parts with the black dashed lines. The temperature threshold is 
2 X 10 " K and the density is split at 10“® cm“®. Each quadrant 
is labeled with the total gas mass inside it. 


a short radial range, the density increases dramatically be¬ 
cause at this density it is able to cool much more efficiently. 
Its cooling time drops below 1 Myr i f the temperature is be - 
low 10^ K. In the terms described in White & Frenk (19911, 
this is the “cooling radius.” 

Moving inside the cooling radius. Fig. [^illustrates the 
transition in gas phase that occurs at n = 10“® cm“®. At 
densities below n < 10“® cm“®, the gas is entirely in a 
warm (> 20000 K), diffuse phase. At higher densities, the 
gas cools, loses pressure and collapses to reach densities of 
0.1-10 cm“®. Images of the cool (T<20000 K), dense gas 
shows that it is in a flattened distribution because it is no 
longer pressure supported, but instead its centripetal accel¬ 
eration supports it. Fig.j^shows that the total gas profile of 
the disc flattens out to though its surface density profile 
is exponential. 

The inner median temperatures evolve with mass. Fig. 
shows that lower mass halos have median temperatures 
of 10 " K in their inner regions. As halo mass increases, the 
median temperature in the inner region initially drops to 10® 
K, but in the highest mass halos, the innermost bins reach 
median temperatures close to 10® K. 

Fig. § is the phase diagram for a high (3 x IO^^Mq) 
mass halo. It shows the multiple gas phases present at high 
densities in the simulations. Three main phases are obvious. 
Below T < 20000 K, Fig. [^ shows two phases of gas, one 
around lO" K and one below 10® K. The branch below 10® 
K has higher metallicity. Higher metallicity gas cools more 
readily, so as halos grow in mass and form more stars, their 
gas is enriched with more chemical elements that facilitate 
cooling to lower temperatures. At these high masses, the 
cooling becomes more effective such that star formation in¬ 
creases. The increase in star formation ejects more energy 
as feedback that creates more hot, dense gas in these simu¬ 
lations. 

The high mass halos have deeper potential wells, so gas 
in the centre of the high mass halos exhibit higher median 


temperatures in Fig. [H While gas of 10® K escapes from 
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lower mass halos, 10^^ M© halos can hold onto this gas. 
Fig. [9] (the phase diagram of a high mass halo) shows that 
there is almost as much mass in the high temperature, high 
density phase as cool gas in the disc. In our simulations, the 
disabled cooling from stellar feedback artificially maintains 
gas in this state. This hot, dense gas is providing the pressure 
support that prevents a major cooling flow. At these stellar 
masses, other observations show that quenching may start 
to set in. So, the stellar feedback is less effective than what is 
happening in the real universe. Thus, many models employ 
strong AGN feedback at these mass scales to decrease star 
formation ( Croton et al.|[2006 Vogelsberger et al.|2014 l. 

One hnal feature to note in Fig. [7| is that lower mass 
halos exhibit a broad scatter in the density profiles of their 
inner regions. The scatter reduces as mass increases. The 
scatter comes from stellar feedback that disrupts the gaseous 
discs. Lower mass halos have shallower potential wells, so it 
is easier for the same amount of stellar feedback energy to 
blow apart the disc. As the galaxies grow in mass, the stellar 
feedback has notably less effect. Fig. shows that enough 
hot gas is created that it becomes the dominant phase in the 
inner regions of the highest mass galaxies in our sample. 


5 DISCUSSION 

5.1 Theoretical Simplification 


Galaxy formation can seem complex and involves many 
physical processes. So, it is comforting when a simple re¬ 
sult comes out of a galaxy formation model. Theorists have 
noticed that if you assume the mass of cold disc gas remains 
constant, it is possible to produce simple models of galaxy 
formation that reproduce many observed propertie s (|Bouche| 
et al. 2010 Dave et al.||2012' Krumholz fc Dekel||2012 1. In 


these models, sometimes referred to as “bathtub models” 
the amount of gas that is accreted onto galaxies is balanced 
by sinks of gas, star formation and outflows to leave a con¬ 
stant gas mass. The reason given for a constant gas mass 
in these models is that the equations governing the models 
include the derivative of gas mass with time. The solution 
to a simple ordinary differential equation is a decaying ex¬ 
ponential, so the time derivative of gas mass should decay 
to 0 ( Dekel fc Mandelker||2014 |. 


5.2 Why Constant Gas Mass? 


The inner region of galaxy halos have mean densities corre¬ 
sponding to dynamical times of 100s Myr. This is less than 
the age of the Universe, so one would expect that some sort 
of steady state could be reached. Specifically, one should 
compare the dynamical and cooling times with the “gas de¬ 
pletion timescale.” 


The “gas depletion timescale” (e.g. Saintonge et al. 


20111 compares the mass of gas observed in discs with their 


star formation rate and outflow rate (the gas sinks). Typical 
observed timescales are ~ 2 Gyr, which is an indication of 
how inefficient star formation is in galaxies: there is a large 
reservoir of cool gas in galaxy discs (~ 10^'^ Mq), but only 
a small fraction of it forms stars (~ 1 Mq yr“^). 

If the cooling and free fall times of halo gas is shorter 
than the gas depletion timescale, then gas will cool into the 


disc and be replenished. Cooling time is a function of tem¬ 
perature and density. As a galaxy grows, its virial temper¬ 
ature increases. At some temperature, the gas will be hot 
enough that its cooling time becomes longer than the gas 
depletion time. At this temperature (~ 10^ K), the hot gas 
will also provide enough pressure to support the hydrostatic 
equilibrium of the galaxy halo. This temperature is beyond 
that which supernova create, but within reach of the tem¬ 
perature gas that outflows from AGN can create. The gas 
depletion timescale then sets the time for “strangulation” of 
the galaxy to occur. 


5.3 Comparing Mgas-M 2 oo with M*-M 2 oo 

Fig. 0 showed that the relation between Mgas and M 200 
has a characteristic shape that is steep, Mgas ~M 2 oo, below 
M 200 = 2 X and then flattens to linear above M 200 > 

2x lO^^M©. In comparison, the M*-Mhaio relation has a sim¬ 
ilarly steep slope, M* oc M 200 for M 200 < 2 x IO^^Mq halos 
at 2 = 0 ( [Kravtsov et al.|2014[ ). 

The slopes of Mgas-M 2 oo and M*-M 2 oo also both evolve 


with redshift. Moster et al. (20131 find that M*-M 2 oo gets 


steeper as time evolves, so that it is steepest at z=0. Fig. 
shows that the same is true for Mgas-M 2 oo. 

Looking above M 200 > 2 x IO^^Mq, there are two cases 


reported for M*-M 2 oo. Moster et al. (20131 report that M*- 
M 200 becomes shallower, while Behroozi et al. ] poT^ find 


Mv.-M;,ai o steepens in the ra nge 2 x 10^^ < M 2 oo/Ma < 
10^^ The 


Moster et al. 


(20131 result indicates that star for¬ 


mation follows the growth of gas mass. The Behroozi result 
points to a mass range in which star formation has enhanced 
efficiency. 


5.4 Observational Verification 


Prochaska et al. {20051 noticed that when you add up all the 


damped Lyman-a systems seen in quasar absorption lines, 
it turned out that the Hi density {Q.hi) in the Universe 
remains constant throughout the history of the Universe. 
We do not have a large enough statistical sample to state 
whether the constant cool disc gas mass corresponds to a 
constant Qhi, but it would not be surprising if there is a 


connection. Interested readers should see Lagos et al. (20111, 


who examines how the cool gas fraction evolves in semi- 
analytic models. 

We have seen in these simulations that the exact treat¬ 
ment of stellar feedback does not greatly affect the steady 
state disc gas mass. Stellar feedback determines the effi¬ 


ciency with which accreting gas turns into stars (a la Hop- 
kins|2013 1. The disc gas mass is set by the conversion of gas 


into stars, the Kennicutt-Schmidt relationship (Kennicutt 
iQ^jBigiel et al.|[2008|. 


White & Frenk ( 1991| note that the exact disc gas mass 


is not important for modeling the galaxy luminosity func¬ 
tion. Gold gas is just one phase gas takes as it cools and turns 
into stars. However, there are some interesting observational 
implications. [Bigiel fc Blitzj ( |2012| ) recently found that iso¬ 
lated galaxies have a universal neutral gas profile that leads 
them to conclude that neutral gas mass only depends on 
their optical radius, r 25 . We are not the only ones seeing a 
correspondence between neutral gas and simulations, jWaii^ 
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et al. (20141 find that the neutral gas in simulations com¬ 


pares well with observations. 


5.5 Why is 3 X 10^°Mq the maximum gas mass? 

To understand what limits the mass of cool gas in discs, It 
is instructional to integrate the mass included in an expo¬ 
nential disc with a typical scale radius. For an exponential 
profile of the form 

E(i?) = Eoe"'^ (3) 


with a scale length, Rd, of 5 kpc, when we integrate out 
to 5 scale lengths (25 kpc), the mass is IO^^Mq when Eq, 
the surface density in the centre, is 100 Mq pc“^, nearly 


the same surface density as typical molecular clouds (Lar- 


1981 McKee & Ostriker 20071. Molecular clouds are 


regions where there is sufficient shielding from external UV 
radiation fields such that molecules will not be dissociated. 
Instead, the molecules radiatively cool gas to low enough 
temperatures to enable collapse to high enough densities for 
fusion to start and stars to form. The maximum disc gas 
mass is the disc mass where star formation becomes effi¬ 
cient. 


5.6 Speculation about quenching star formation 

Since our simulated sample is limited at the high mass end 
as the simulations overproduce stars in the galactic cen¬ 
tres, we are forced to speculate about how star formation is 
quenched. Extrapolating our results indicates that increas¬ 
ing virial temperatures possibly combined with other ener¬ 
getic input (AGN?) raise cooling times sufficiently to halt 
gas accretion onto galaxies and quench star formation. We 
see a slight indication of decreasing cool gas mass in our 
highest mass simulations, but those simulations do not re¬ 
duce their gas mass quickly enough to correspond to galaxy 
quenching. 


6 CONCLUSIONS 

We study the gas masses in a large set of cosmological hy¬ 
drodynamic simulations that well reproduce both the halo 
M*-M 2 oo relation and the M*-Mhi relation in observations. 
We find that disc gas masses remain remarkably constant 
over many Gyrs following their last major merger. The cool 
disc masses vary less than 50% in all galaxies such that most 
galaxy cool gas masses vary by less than unity over the Hub¬ 
ble time. Since cool gas masses stay constant for long periods 
of time, the steady state of galaxies can be modeled using 
equilibrium or “bathtub” models. 

The most important factor in how cool gas mass evolves 
is how fast the halo spins. Low spin galaxies decrease their 
cool disc gas mass while high spin galaxies can add cool 
gas. The difference corresponds to the efficiency of star for¬ 
mation. Low spin galaxies reach higher gas surface densities 
that lead to more efficient star formation. Consequently, low 
spin galaxies have lower disc gas fractions than high spin 
galaxies. 

Over a wide range of galaxy masses, the constant cool 
disc gas mass follows a tight correlation with the halo virial 


mass. The disc gas mass-halo mass relation has a break at 
M 200 = 2 X IO^'^Mq. At masses below the break, Mgas ~ 
M 200 ; above the break, Mgas ~ M^oo- The low mass rela¬ 
tion steepens with time, likely in response to changes in the 
ionization background. Observations show a similar break 
in the Mgas-M* relation, and a comparison shows that the 
simulations break at a similar stellar mass to the observa¬ 
tions. 

The tight relation between halo mass and cool gas mass 
indicates that there is some sort of universal physical mech¬ 
anism that determines the cool gas mass. To study the 
mechanism that sets the cool gas mass, we examine the 
one-dimensional total gas density profiles of galaxies. The 
profiles show three characteristic features across our entire 
sample. In the outskirts, the gas halos follow an isothermal 
profile, p oc r~^. When the density reaches 10“^ cm“®, the 
density steeply increases to 0.1 cm“^ due to the shortened 
cooling time in dense gas. Inside that region, the ID spher¬ 
ical density profile of the disc is a flatter p oc r~^ power 
law, though it shows large scatter in the lowest mass halos. 
The scattered density profile results from stellar feedback 
disrupting the disc in the lowest mass halos. 

The constant gas masses and near universal density pro¬ 
files lead us to conclude that gas in galaxies reaches a quasi¬ 
static equilibrium configuration. In this configuration, star 
formation and outflows result from any accretion onto the 
disc. Star formation rates and outflows are relatively small 
due to the large amount of gas mass contained in the diffuse 
gas halo. Thus, stellar feedback must be effective at pro¬ 
viding pressure support to prevent mass cooling flows from 
hot halo gas onto the disc for galaxies with lower mass than 
the Milky Way. We leave studying how the high entropy 
(low density, high temperature) gas that stellar feedback 
produces interacts with the halo for future work. 

In a small experiment that explored how the constant 
gas mass is affected when various feedbacks are employed, 
we find that feedback does not have a large effect on the 
constant gas mass. While the gas disc masses are relatively 
constant, the mass of stars in each simulation vary greatly 
depending on the specifics of the star formation and feedback 
recipe. How much cool gas exists in the disc should depend 
sensitively on the difficult problem of correctly modeling the 
stark boundary between cold, dense gas in the galaxy disc 
and the hot, diffuse gas in the halo surrounding galaxies. 
Thus, cool gas mass should be a good metric to compare 
how various hydrodynamic methods simulate galaxies. 
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